Binary responses of corals to sediment stress
Lillian J. Tuttle
9/24/2020
Here I will explore how DEPOSITED AND SUSPENDED sediment affects the PRESENCE/ABSENCE OF NON-ADVERSE AND ADVERSE EFFECTS of tropical, scleractinean coral of all life-history stages. The database to be used is comprised of data extracted from studies deemed relevant during the systematic literature review process.
The Dataset
Deposited Sediment
Specifically, we will explore the results of 1644 datapoints from 45 studies conducted in 3 oceans on 113 species within 53 genera (shown below only the top 10 species and genera, by number of datapoints).
## n
## 1 1644## Ref n
## 1 DS01 24
## 2 DS02 16
## 3 DS03 12
## 4 DS04 70
## 5 DS05 10
## 6 DS06 6
## 7 DS07 13
## 8 DS08a 162
## 9 DS08b 113
## 10 DS10 16
## 11 DS11 4
## 12 DS12 7
## 13 DS13 126
## 14 DS15 54
## 15 DS16 46
## 16 DS17 24
## 17 DS18 16
## 18 DS19 62
## 19 DS20 8
## 20 DS21 4
## 21 DS23 51
## 22 DS24 4
## 23 DS25 2
## 24 DS26 6
## 25 DS27 39
## 26 DS28 4
## 27 DS29 40
## 28 DS30 36
## 29 DS31 168
## 30 DS32 16
## 31 DS33 8
## 32 DS34 157
## 33 DS36 2
## 34 DS37 4
## 35 DS38 24
## 36 DS40 20
## 37 DS42 2
## 38 DS43 12
## 39 DS45 24
## 40 DS46 44
## 41 DS48 20
## 42 DS49 6
## 43 DS68 48
## 44 DS69 2
## 45 DS71 112## Ref Ref_name
## 1 DS01 Babcock and Davies (1991)
## 25 DS02 Babcock and Smith (2000)
## 41 DS03 Bessell-Browne et al. (2017a)
## 53 DS04 Duckworth et al. (2017)
## 123 DS05 Hodel (2007)
## 133 DS06 Hodgson (1990a)
## 139 DS07 Hodgson (1990b)
## 152 DS08a Hodgson (1989) PhD Chapter IV
## 314 DS08b Hodgson (1989) PhD Chapter IV
## 427 DS10 Junjie et al. (2014)
## 443 DS11 Lirman et al. (2008)
## 447 DS12 Loiola et al. (2013)
## 454 DS13 Moeller et al. (2017)
## 580 DS15 Perez et al. (2014)
## 634 DS16 Philipp and Fabricius (2003)
## 680 DS17 Piniak (2007)
## 704 DS18 Piniak and Brown (2008)
## 720 DS19 Ricardo et al. (2017)
## 782 DS20 Riegl and Branch (1995)
## 790 DS21 Rogers (1979)
## 794 DS23 Selim (2007)
## 845 DS24 Sheridan et al. (2014)
## 849 DS25 Shore-Maggio et al. (2018)
## 851 DS26 Sofonia (2006) PhD Chapter 3
## 857 DS27 Sofonia (2006) PhD Chapter 4
## 896 DS28 Sofonia and Anthony (2008)
## 900 DS29 Stafford-Smith (1990) PhD Chapter 4
## 940 DS30 Stafford-Smith (1992)
## 976 DS31 Stafford-Smith and Ormond (1992)
## 1144 DS32 Stewart et al. (2006)
## 1160 DS33 Vargas-Angel et al. (2006)
## 1168 DS34 Weber et al. (2006)
## 1325 DS36 Zill et al. (2017)
## 1327 DS37 Fabricius et al. (2003)
## 1331 DS38 Goh and Lee (2008)
## 1355 DS40 Abdel-Salam (1989) PhD Chapter 3
## 1375 DS42 Gowan et al. (2014)
## 1377 DS43 Peters and Pilson (1985)
## 1389 DS45 Rogers (1983)
## 1413 DS46 Stafford-Smith (1993)
## 1457 DS48 Bessell-Browne et al. (2017) Lab component
## 1477 DS49 Coffroth (1985)
## 1483 DS68 Flores et al. (2012)
## 1531 DS69 Gil et al. (2016)
## 1533 DS71 HDR EOC and CSA Ocean Services, Inc. (2014)## Ref_name n
## 1 Abdel-Salam (1989) PhD Chapter 3 20
## 2 Babcock and Davies (1991) 24
## 3 Babcock and Smith (2000) 16
## 4 Bessell-Browne et al. (2017) Lab component 20
## 5 Bessell-Browne et al. (2017a) 12
## 6 Coffroth (1985) 6
## 7 Duckworth et al. (2017) 70
## 8 Fabricius et al. (2003) 4
## 9 Flores et al. (2012) 48
## 10 Gil et al. (2016) 2
## 11 Goh and Lee (2008) 24
## 12 Gowan et al. (2014) 2
## 13 HDR EOC and CSA Ocean Services, Inc. (2014) 112
## 14 Hodel (2007) 10
## 15 Hodgson (1989) PhD Chapter IV 275
## 16 Hodgson (1990a) 6
## 17 Hodgson (1990b) 13
## 18 Junjie et al. (2014) 16
## 19 Lirman et al. (2008) 4
## 20 Loiola et al. (2013) 7
## 21 Moeller et al. (2017) 126
## 22 Perez et al. (2014) 54
## 23 Peters and Pilson (1985) 12
## 24 Philipp and Fabricius (2003) 46
## 25 Piniak (2007) 24
## 26 Piniak and Brown (2008) 16
## 27 Ricardo et al. (2017) 62
## 28 Riegl and Branch (1995) 8
## 29 Rogers (1979) 4
## 30 Rogers (1983) 24
## 31 Selim (2007) 51
## 32 Sheridan et al. (2014) 4
## 33 Shore-Maggio et al. (2018) 2
## 34 Sofonia (2006) PhD Chapter 3 6
## 35 Sofonia (2006) PhD Chapter 4 39
## 36 Sofonia and Anthony (2008) 4
## 37 Stafford-Smith (1990) PhD Chapter 4 40
## 38 Stafford-Smith (1992) 36
## 39 Stafford-Smith (1993) 44
## 40 Stafford-Smith and Ormond (1992) 168
## 41 Stewart et al. (2006) 16
## 42 Vargas-Angel et al. (2006) 8
## 43 Weber et al. (2006) 157
## 44 Zill et al. (2017) 2## Ocean n
## 1 Pacific 1446
## 2 Atlantic 95
## 3 Indian 63
## 4 Indo-Pacific 40## Selecting by n## Gsp n
## 1 Montipora_peltiformis 167
## 2 Acropora_millepora 140
## 3 Pocillopora_damicornis 121
## 4 Leptastrea_purpurea 99
## 5 Porites_lobata/lutea 86
## 6 Acropora_hyacinthus 50
## 7 Leptoria_phrygia 46
## 8 Porites_cylindrica 43
## 9 Montipora_aequituberculata 30
## 10 Porites_rus 30## Selecting by n## Updated_Genus n
## 1 Montipora 319
## 2 Acropora 258
## 3 Porites 205
## 4 Pocillopora 133
## 5 Leptastrea 99
## 6 Turbinaria 47
## 7 Leptoria 46
## 8 Pavona 37
## 9 Echinopora 31
## 10 Oxypora 31Suspended Sediment
Specifically, we will explore the results of 661 datapoints from 41 studies conducted in 3 oceans on 29 species within 18 genera (shown below only the top 10 species and genera, by number of datapoints).
## n
## 1 661## Ref n
## 1 SS06 150
## 2 SS08 48
## 3 SS14c 43
## 4 SS20d 37
## 5 SS16a 32
## 6 SS22c 28
## 7 SS17b 26
## 8 SS17a 24
## 9 SS21a 23
## 10 SS15 20
## 11 SS16c 16
## 12 SS27 16
## 13 SS20c 15
## 14 SS01 12
## 15 SS03a 12
## 16 SS05 12
## 17 SS11a 12
## 18 SS12a 12
## 19 SS22a 10
## 20 SS04 9
## 21 SS12b 8
## 22 SS16b 8
## 23 SS20a 8
## 24 SS22b 8
## 25 SS19a 6
## 26 SS20b 6
## 27 SS21b 6
## 28 SS13a 5
## 29 SS13b 5
## 30 SS13c 5
## 31 SS13d 5
## 32 SS14a 5
## 33 SS19b 4
## 34 SS24a 4
## 35 SS24b 4
## 36 SS07 3
## 37 SS11b 3
## 38 SS11c 3
## 39 SS14b 3
## 40 SS25 3
## 41 SS28 2## Ref_name n
## 1 Browne et al. 2015 150
## 2 Ricardo et al. 2016 66
## 3 Kendall et al. 1985 56
## 4 Humphrey et al. 2008 51
## 5 Liu et al. 2015 50
## 6 Flores et al. 2012 48
## 7 Rice 1984 46
## 8 Ricardo et al. 2018 29
## 9 Humanes et al. 2017a 20
## 10 Humanes et al. 2017b 20
## 11 Jokiel et al. 2014 20
## 12 Gilmour 1999 18
## 13 Anthony et al. 2007 16
## 14 Anthony 1999 12
## 15 Anthony and Fabricius 2000 12
## 16 Browne et al. 2014 12
## 17 Ricardo et al. 2015 10
## 18 Bessell-Browne et al. 2017 9
## 19 Te 1992 8
## 20 Erftemeijer et al. 2012 3
## 21 Te 2001 Chapter 6 3
## 22 Dallmeyer et al. 1982 2## Ocean n
## 1 Pacific 392
## 2 Indian 153
## 3 Atlantic 104
## 4 Indo-Pacific (e.g. Singapore) 12## Selecting by n## Gsp n
## 1 Acropora_millepora 131
## 2 Acropora_tenuis 77
## 3 Acropora_cervicornis 56
## 4 Merulina_ampliata 54
## 5 Pachyseris_speciosa 54
## 6 Platygyra_sinensis 54
## 7 Acropora_muricata 50
## 8 Montipora_aequituberculata 30
## 9 Acropora_digitifera 18
## 10 Acropora_intermedia 16## Selecting by n## Updated_Genus n
## 1 Acropora 348
## 2 Merulina 54
## 3 Pachyseris 54
## 4 Platygyra 54
## 5 Montipora 46
## 6 Porites 25
## 7 Pocillopora 17
## 8 Goniastrea 12
## 9 Cladocora 8
## 10 Manicina 8
## 11 Solenastrea 8Exploratory plots
First let’s explore all data from all species for which there is binary data about the presence of ‘non-adverse effects’, ‘adverse effects’, and ‘any mortality’ as a result of exposure to deposited and suspended sediment.
DEFINITIONS:
- ‘Non-adverse effects’ are defined as hydrostatic inflation, tentacular activity, increased mucus production, and/or congealed sediment.
- ‘Adverse effects’ are defined as a shift from auto- to hetero-trophy, reduced photosynthetic efficiency, reduced P/R ratio, reduced growth rate, bleaching (or reduced chlorophyll A content), algal overgrowth, smothering, and any degree of tissue/colony mortality. Here we consider only adult corals.
- ‘Any mortality’ is defined as the presence of small, large, or total necroses of coral’s surface area.
Threshold Analyses with Binary Data
Now let’s calculate thresholds based on the binary data explored above.
DEFINITIONS:
- ‘LOAEL’, or the ‘lowest observed adverse effect level’ is defined as the lowest dose at which there was an observed adverse effect.
- ‘NOAEL’, or the ‘no observed adverse effect level’ is defined as the highest dose at which there was NOT an observed adverse effect.
- ‘Adverse effect’ is defined here as any response of a coral individual, colony, or experimental treatment group that may negatively affect the coral’s fitness and/or survival. These adverse effects may include sublethal physiological changes (e.g., significantly reduced growth or photosynthetic rates when compared to coral in ambient/control conditions), bleaching/tissue loss, and mortality. This definition of adverse effect is independent of its magnitude, so while the affect may have potentially reduced fitness, its magnitude may be sufficiently small that the fitness effect is not measurable.
DEPOSITED SEDIMENT
Adverse effect
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 1## NOAEL
## 1 1
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 0.5## NOAEL
## 1 0.5
Any mortality
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 4.9## NOAEL
## 1 4.4
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 0.917## NOAEL
## 1 0.917
SUSPENDED SEDIMENT
Adverse effect
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 3.22## NOAEL
## 1 3.22
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 0.04166667## NOAEL
## 1 0.0208333
Any mortality
NOAEL/LOAEL defined by exposure concentration
## LOAEL
## 1 3.22## NOAEL
## 1 3.22
NOAEL/LOAEL defined by exposure duration
## LOAEL
## 1 0.5## NOAEL
## 1 0.5
Logistic Regression Model Fitting
DEPOSITED SEDIMENT
Adverse effect
## n
## 1 943## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Data: adverseDS3
## Models:
## modDS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS6: (1 | Gsp)
## modDS9: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS9: (1 | Ref)
## modDS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## modDS15: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS15: (1 | Ref_name/Ref)
## modDS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS18: (1 | Gsp) + (1 | Ref_name)
## modDS21: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS21: (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS6 5 1129.5 1153.8 -559.75 1119.5
## modDS9 5 1234.8 1259.1 -612.42 1224.8 0.0000 0 1.0000
## modDS12 6 1078.3 1107.3 -533.13 1066.3 158.5888 1 <2e-16 ***
## modDS15 6 1234.9 1264.0 -611.45 1222.9 0.0000 0 1.0000
## modDS18 6 1084.3 1113.4 -536.16 1072.3 150.5848 0 <2e-16 ***
## modDS21 7 1080.1 1114.0 -533.04 1066.1 6.2525 1 0.0124 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: adverseDS3
## Models:
## modDS10: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDS10: Ref)
## modDS11: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days +
## modDS11: (1 | Gsp) + (1 | Ref)
## modDS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modDS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS10 4 1106.5 1125.9 -549.26 1098.5
## modDS11 5 1085.9 1110.1 -537.93 1075.9 22.6545 1 1.939e-06 ***
## modDS12 6 1078.3 1107.3 -533.13 1066.3 9.5959 1 0.00195 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: adverseDS3
## Models:
## modDS4: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp)
## modDS7: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Ref)
## modDS10: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDS10: Ref)
## modDS13: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modDS16: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDS16: Ref_name)
## modDS19: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDS19: Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDS4 3 1133.9 1148.4 -563.93 1127.9
## modDS7 3 1244.0 1258.6 -619.01 1238.0 0.00 0 1.0000
## modDS10 4 1106.5 1125.9 -549.26 1098.5 139.50 1 <2e-16 ***
## modDS13 4 1242.3 1261.7 -617.13 1234.3 0.00 0 1.0000
## modDS16 4 1104.2 1123.6 -548.11 1096.2 138.03 0 <2e-16 ***
## modDS19 5 1106.2 1130.5 -548.11 1096.2 0.00 1 0.9991
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days +
## (1 | Gsp) + (1 | Ref)
## Data: adverseDS3
##
## AIC BIC logLik deviance df.resid
## 1085.9 1110.1 -537.9 1075.9 938
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0176 -0.6502 -0.1038 0.5777 3.4105
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 5.151 2.270
## Ref (Intercept) 4.313 2.077
## Number of obs: 943, groups: Gsp, 91; Ref, 38
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.233170 0.544485 -2.265 0.023523 *
## Sed_level_stand_mg 0.004003 0.001202 3.330 0.000867 ***
## Sed_exposure_days 0.019949 0.004590 4.346 1.38e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_effect_adverse ~ log10_conc + Sed_exposure_days + (1 |
## Gsp) + (1 | Ref)
## Data: adverseDS3
##
## AIC BIC logLik deviance df.resid
## 1067.9 1092.2 -529.0 1057.9 938
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8148 -0.6181 -0.0858 0.5431 3.4637
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 5.470 2.339
## Ref (Intercept) 5.289 2.300
## Number of obs: 943, groups: Gsp, 91; Ref, 38
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.019371 0.711298 -4.245 2.19e-05 ***
## log10_conc 1.236537 0.239249 5.168 2.36e-07 ***
## Sed_exposure_days 0.021173 0.004644 4.559 5.13e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.804878##
## p 0 1
## 0 412 116
## 1 68 347## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.8877## R2m R2c
## theoretical 0.05186506 0.7779635
## delta 0.04946068 0.7418984
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 2.3387
## Ref (Intercept) 2.2998## [1] 9.971703
## Est LL UL
## (Intercept) 0.04883195 0.01211253 0.1968671
## log10_conc 3.44366632 2.15460790 5.5039424
## Sed_exposure_days 1.02139871 1.01214419 1.0307379## Est LL UL
## (Intercept) 9.563777e-04 3.859093e-05 0.02370138
## log10_conc 1.723998e+01 5.856104e+00 50.75334043
## Sed_exposure_days 1.049961e+00 1.028184e+00 1.07219781For every 10-fold increase in exposure concentration of deposited sediment, the odds of an adverse effect increase by 3.4 times (95% CI 2.2, 5.5, GLMM z = 5.168, p < 0.0001), while holding exposure duration constant and accounting for variability among studies and species. For every 10-fold increase in exposure duration of deposited sediment, the odds of an adverse effect increase by 1.05 times (95% CI 1.03, 1.07, GLMM z = 4.559, p < 0.0001), while holding exposure concentration constant and accounting for variability among studies and species.
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS11_log <- predict(modDS11_log, newdata=adverseDS3, type="response")
adverseDS3 <- cbind(adverseDS3, pred_modDS11_log)
adverseDS_plot <- ggplot(data = adverseDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_adverse)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Adverse effects due to sediment exposure?") +
scale_x_log10(limits=c(0.01,max(adverseDS3$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS11_log), inherit.aes=FALSE) +
theme_bw()
adverseDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(adverseDS3, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
adverseDS3$log10_conc <- j
predict(modDS11_log, newdata = adverseDS3, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(0,1,2,3)], mean)
##
# get the means with lower and upper quartiles
plotdat1_2 <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat1_2 <- as.data.frame(cbind(plotdat1_2, jvalues))
plotdat1_2 <- plotdat1_2 %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat1_2) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat1_2)
#approx(plotdat1_2$Sed_conc, plotdat1_2$PredictedProbability, xout=1)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat1_2$Sed_conc, plotdat1_2$PredictedProbability, xout=x))
unlist(sed2)## x y x y x y
## 1.0000000 0.1922209 5.0000000 0.3039429 10.0000000 0.3593714
## x y x y x y
## 50.0000000 0.4971127 100.0000000 0.5567901 500.0000000 0.6854517
## x y
## 1000.0000000 0.7335196# plot average marginal predicted probabilities
ggplot(plotdat1_2, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat1_2, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat1_2, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
adverseDS_plot2 <- ggplot() +
geom_point(data = adverseDS3, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_adverse)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of adverse\neffect due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.01,max(adverseDS3$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000,loael),
label=c("0.01","0.1","1","10","100","1000",round(loael,digits=1))) +
geom_ribbon(data = plotdat1_2, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat1_2, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat1_2, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loael, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
adverseDS_plot2
Any mortality
## n
## 1 719## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Data: anymortDS2
## Models:
## modDSb6: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb6: (1 | Gsp)
## modDSb9: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb9: (1 | Ref)
## modDSb12: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## modDSb15: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb15: (1 | Ref_name/Ref)
## modDSb18: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb18: (1 | Gsp) + (1 | Ref_name)
## modDSb21: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb21: (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb6 5 758.73 781.62 -374.36 748.73
## modDSb9 5 840.41 863.30 -415.20 830.41 0.0000 0 1.000000
## modDSb12 6 733.43 760.90 -360.72 721.43 108.9770 1 < 2.2e-16 ***
## modDSb15 6 839.84 867.31 -413.92 827.84 0.0000 0 1.000000
## modDSb18 6 742.10 769.56 -365.05 730.10 97.7451 0 < 2.2e-16 ***
## modDSb21 7 735.05 767.10 -360.53 721.05 9.0432 1 0.002637 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: anymortDS2
## Models:
## modDSb10: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDSb10: Ref)
## modDSb11: Binary_effect_mortality ~ Sed_level_stand_mg + Sed_exposure_days +
## modDSb11: (1 | Gsp) + (1 | Ref)
## modDSb12: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modDSb12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb10 4 768.46 786.77 -380.23 760.46
## modDSb11 5 738.10 760.99 -364.05 728.10 32.3569 1 1.283e-08 ***
## modDSb12 6 733.43 760.90 -360.72 721.43 6.6698 1 0.009806 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: anymortDS2
## Models:
## modDSb4: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp)
## modDSb7: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Ref)
## modDSb10: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDSb10: Ref)
## modDSb13: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modDSb16: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDSb16: Ref_name)
## modDSb19: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modDSb19: Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modDSb4 3 790.99 804.73 -392.50 784.99
## modDSb7 3 855.44 869.17 -424.72 849.44 0.000 0 1.0000
## modDSb10 4 768.46 786.77 -380.23 760.46 88.981 1 <2e-16 ***
## modDSb13 4 852.46 870.77 -422.23 844.46 0.000 0 1.0000
## modDSb16 4 765.01 783.32 -378.51 757.01 87.445 0 <2e-16 ***
## modDSb19 5 767.01 789.90 -378.51 757.01 0.000 1 0.9997
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_effect_mortality ~ log10_conc + Sed_exposure_days + (1 |
## Gsp) + (1 | Ref)
## Data: anymortDS2
##
## AIC BIC logLik deviance df.resid
## 721.3 744.2 -355.6 711.3 714
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3408 -0.4784 -0.0902 0.4327 3.2066
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 7.998 2.828
## Ref (Intercept) 5.484 2.342
## Number of obs: 719, groups: Gsp, 81; Ref, 30
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.749716 1.251060 -5.395 6.84e-08 ***
## log10_conc 2.139224 0.444329 4.814 1.48e-06 ***
## Sed_exposure_days 0.032171 0.006279 5.123 3.00e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.8484006##
## p 0 1
## 0 363 56
## 1 53 247## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9317## R2m R2c
## theoretical 0.09992789 0.8234467
## delta 0.09445315 0.7783326
## $Gsp
##
## $Ref
## Groups Name Std.Dev.
## Gsp (Intercept) 2.8280
## Ref (Intercept) 2.3419## [1] 10.40077
## Est LL UL
## (Intercept) 0.001171212 0.0001008583 0.01360065
## log10_conc 8.492841038 3.5549391862 20.28961541
## Sed_exposure_days 1.032693755 1.0200618185 1.04548212## Est LL UL
## (Intercept) 1.779443e-07 6.283564e-10 5.039204e-05
## log10_conc 1.377919e+02 1.854975e+01 1.023550e+03
## Sed_exposure_days 1.076888e+00 1.046799e+00 1.107843e+00For every 10-fold increase in exposure concentration of deposited sediment, the odds of any tissue mortality increase by 8.5 times (95% CI 3.6, 20.3; GLMM z = 4.814, p < 0.0001), while holding exposure duration constant and accounting for variability among studies and species. For every 10-fold increase in exposure duration of deposited sediment, the odds of any tissue mortality increase by 1.08 times (95% CI 1.05, 1.11; GLMM z = 5.123, p < 0.0001), while holding exposure concentration constant and accounting for variability among studies and species.
Prediction figure
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSb11_log <- predict(modDSb11_log, newdata=anymortDS2, type="response")
anymortDS3 <- cbind(anymortDS2, pred_modDSb11_log)
anymortDS_plot <- ggplot(data = anymortDS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_mortality,
color = Ref)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Any mortality due to sediment exposure?",
color = "Study") +
scale_x_log10(limits=c(0.01,1000),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSb11_log), inherit.aes=FALSE) +
theme_bw()
anymortDS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
jvalues <- with(anymortDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
anymortDS2$log10_conc <- j
predict(modDSb11_log, newdata = anymortDS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
##
# get the means with lower and upper quartiles
plotdat2_2 <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat2_2 <- as.data.frame(cbind(plotdat2_2, jvalues))
plotdat2_2 <- plotdat2_2 %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat2_2) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat2_2)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat2_2$Sed_conc, plotdat2_2$PredictedProbability, xout=x))
print(unlist(sed2), digits = 3)## x y x y x y x y
## 1.00e+00 7.96e-02 5.00e+00 1.92e-01 1.00e+01 2.58e-01 5.00e+01 4.40e-01
## x y x y x y
## 1.00e+02 5.21e-01 5.00e+02 6.74e-01 1.00e+03 7.25e-01# plot average marginal predicted probabilities
ggplot(plotdat2_2, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat2_2, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat2_2, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
anymortDS_plot2 <- ggplot() +
geom_point(data = anymortDS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_mortality)) +
labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"),
y = "Predicted probability of any tissue\nmortality due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.01,1000), breaks=c(0.01,0.1,1,10,100,1000,loael2),
label=c("0.01","0.1","1","10","100","1000",round(loael2,digits=1))) +
geom_ribbon(data = plotdat2_2, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat2_2, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat2_2, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loael2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
anymortDS_plot2
Suspended Sediment
Adverse Effect
## n
## 1 423## Ref n
## 1 SS01 8
## 2 SS03a 8
## 3 SS04 6
## 4 SS05 9
## 5 SS06 108
## 6 SS07 2
## 7 SS08 30
## 8 SS11a 8
## 9 SS11b 2
## 10 SS11c 2
## 11 SS12a 9
## 12 SS12b 6
## 13 SS13a 4
## 14 SS13b 4
## 15 SS13c 4
## 16 SS13d 4
## 17 SS14a 4
## 18 SS14b 2
## 19 SS14c 38
## 20 SS17a 20
## 21 SS17b 22
## 22 SS19a 3
## 23 SS19b 2
## 24 SS20a 7
## 25 SS20b 4
## 26 SS20c 13
## 27 SS20d 32
## 28 SS21a 20
## 29 SS21b 3
## 30 SS22a 5
## 31 SS22b 4
## 32 SS22c 14
## 33 SS24a 3
## 34 SS24b 3
## 35 SS25 3
## 36 SS27 6
## 37 SS28 1## Ref_name n
## 1 Anthony 1999 8
## 2 Anthony and Fabricius 2000 8
## 3 Anthony et al. 2007 6
## 4 Bessell-Browne et al. 2017 6
## 5 Browne et al. 2014 9
## 6 Browne et al. 2015 108
## 7 Dallmeyer et al. 1982 1
## 8 Erftemeijer et al. 2012 2
## 9 Flores et al. 2012 30
## 10 Gilmour 1999 12
## 11 Humanes et al. 2017a 15
## 12 Humanes et al. 2017b 16
## 13 Humphrey et al. 2008 44
## 14 Liu et al. 2015 42
## 15 Ricardo et al. 2015 5
## 16 Ricardo et al. 2016 56
## 17 Ricardo et al. 2018 23
## 18 Rice 1984 23
## 19 Te 1992 6
## 20 Te 2001 Chapter 6 3## Data: adverseSS3
## Models:
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS9: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS9: (1 | Ref)
## modSS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS15: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS15: (1 | Ref_name/Ref)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS21: (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 358.90 379.14 -174.45 348.90
## modSS9 5 390.53 410.77 -190.27 380.53 0.000 0 1.0000
## modSS12 6 357.26 381.55 -172.63 345.26 35.270 1 2.871e-09 ***
## modSS15 6 378.86 403.14 -183.43 366.86 0.000 0 1.0000
## modSS18 6 357.16 381.44 -172.58 345.16 21.701 0 < 2.2e-16 ***
## modSS21 7 357.29 385.62 -171.64 343.29 1.870 1 0.1715
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: adverseSS3
## Models:
## modSS4: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp)
## modSS5: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS5: (1 | Gsp)
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS4 3 393.82 405.97 -193.91 387.82
## modSS5 4 362.91 379.10 -177.46 354.91 32.9146 1 9.63e-09 ***
## modSS6 5 358.90 379.14 -174.45 348.90 6.0078 1 0.01424 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: adverseSS3
## Models:
## modSS10: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSS10: Ref)
## modSS11: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS11: (1 | Gsp) + (1 | Ref)
## modSS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS10 4 385.49 401.68 -188.74 377.49
## modSS11 5 362.14 382.38 -176.07 352.14 25.3477 1 4.787e-07 ***
## modSS12 6 357.26 381.55 -172.63 345.26 6.8734 1 0.008749 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: adverseSS3
## Models:
## modSS16: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSS16: Ref_name)
## modSS17: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS17: (1 | Gsp) + (1 | Ref_name)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS16 4 382.83 399.02 -187.42 374.83
## modSS17 5 361.73 381.96 -175.86 351.73 23.1040 1 1.535e-06 ***
## modSS18 6 357.16 381.44 -172.58 345.16 6.5686 1 0.01038 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: adverseSS3
## Models:
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 358.90 379.14 -174.45 348.90
## modSS12 6 357.26 381.55 -172.63 345.26 3.6365 1 0.05652 .
## modSS18 6 357.16 381.44 -172.58 345.16 0.1061 0 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: adverseSS3
## Models:
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 358.90 379.14 -174.45 348.90
## modSS18 6 357.16 381.44 -172.58 345.16 3.7427 1 0.05304 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_effect_adverse ~ log10_conc * Sed_exposure_days + (1 |
## Gsp) + (1 | Ref_name)
## Data: adverseSS3
##
## AIC BIC logLik deviance df.resid
## 337.0 361.3 -162.5 325.0 417
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7039 -0.3983 -0.1620 -0.0356 9.8721
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 9.4806 3.0791
## Ref_name (Intercept) 0.3828 0.6187
## Number of obs: 423, groups: Gsp, 26; Ref_name, 20
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.698662 1.458821 -5.277 1.31e-07 ***
## log10_conc 1.996090 0.425708 4.689 2.75e-06 ***
## Sed_exposure_days 0.047097 0.018551 2.539 0.0111 *
## log10_conc:Sed_exposure_days 0.006267 0.011124 0.563 0.5732
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.8534279##
## p 0 1
## 0 304 43
## 1 19 57## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9084## R2m R2c
## theoretical 0.1372417 0.7842088
## delta 0.1215667 0.6946405
## $Gsp
##
## $Ref_name
## Groups Name Std.Dev.
## Gsp (Intercept) 3.07906
## Ref_name (Intercept) 0.61869## [1] 1.856486
## Est LL UL
## (Intercept) 0.0004534332 2.598605e-05 0.007912004
## log10_conc 7.3602211982 3.195370e+00 16.953547607
## Sed_exposure_days 1.0482241823 1.010796e+00 1.087038738
## log10_conc:Sed_exposure_days 1.0062864184 9.845841e-01 1.028467094## Est LL UL
## (Intercept) 2.001417e-08 2.767255e-11 1.447524e-05
## log10_conc 9.910373e+01 1.451118e+01 6.768261e+02
## Sed_exposure_days 1.114545e+00 1.025033e+00 1.211873e+00
## log10_conc:Sed_exposure_days 1.014534e+00 9.648594e-01 1.066767e+00For every 10-fold increase in exposure concentration of suspended sediment, the odds of experiencing an adverse effect increase by 7.4 times (95% CI 3.2, 17.0; GLMM z = 4.689, p < 0.0001), while holding constant exposure duration and the interaction between concentration and duration, and while accounting for variability among studies and species. For every 10-fold increase in exposure duration of suspended sediment, the odds of experiencing an adverse effect increase by 1.1 times (95% CI 1.0, 1.2; GLMM z = 2.539, p = 0.011), also while holding constant exposure concentration and the interaction between concentration and duration, and while accounting for variability among studies and species.
Prediction figures
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS18_log <- predict(modSS18_log, newdata=adverseSS3, type="response")
adverseSS3 <- cbind(adverseSS3, pred_modSS18_log)
adverseSS_plot <- ggplot(data = adverseSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_adverse,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Adverse effect due to sediment exposure",
color = "Study") +
scale_x_log10(limits=c(0.1,max(adverseSS3$Sed_level_stand_mg)),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS18_log), inherit.aes=FALSE) +
theme_bw()
adverseSS_plot
That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By sediment exposure concentration
jvalues <- with(adverseSS3, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
adverseSS3$log10_conc <- j
predict(modSS18_log, newdata = adverseSS3, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat3_2 <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat3_2 <- as.data.frame(cbind(plotdat3_2, jvalues))
plotdat3_2 <- plotdat3_2 %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat3_2) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat3_2)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat3_2$Sed_conc, plotdat3_2$PredictedProbability, xout=x))
print(round(unlist(sed2), digits=3))## x y x y x y x y
## 1.000 0.022 5.000 0.068 10.000 0.100 50.000 0.212
## x y x y x y
## 100.000 0.280 500.000 0.483 1000.000 0.580# plot average marginal predicted probabilities
ggplot(plotdat3_2, aes(x = log10_conc)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat3_2, aes(x = log10_conc)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat3_2, aes(x = Sed_conc, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
adverseSS_plot2 <- ggplot() +
geom_point(data = adverseSS3, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_adverse)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of adverse\neffect due to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,max(adverseSS3$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000,loaelSS),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS,digits=1))) +
geom_ribbon(data = plotdat3_2, aes(x = Sed_conc, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat3_2, aes(x = Sed_conc, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat3_2, aes(x = Sed_conc, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
adverseSS_plot2
By sediment exposure duration
jvalues <- with(adverseSS3, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
adverseSS3$Sed_exposure_days <- j
predict(modSS18_log, newdata = adverseSS3, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat3_2b <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat3_2b <- as.data.frame(cbind(plotdat3_2b, jvalues))
# better names and show the first few rows
colnames(plotdat3_2b) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat3_2b)
# plot average marginal predicted probabilities
ggplot(plotdat3_2b, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat3_2b, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
adverseSS_plot3 <- ggplot() +
geom_point(data = adverseSS3, mapping = aes(
x = Sed_exposure_days,
y = Binary_effect_adverse,
color = Ref)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of adverse\neffect due to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(adverseSS3$Sed_exposure_days))) +
geom_ribbon(data = plotdat3_2b, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat3_2b, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat3_2b, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
adverseSS_plot3
Any mortality
## n
## 1 261## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Data: anymortSS2
## Models:
## modSS6: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS6: (1 | Gsp)
## modSS9: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS9: (1 | Ref)
## modSS12: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS12: (1 | Gsp) + (1 | Ref)
## modSS15: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS15: (1 | Ref_name/Ref)
## modSS18: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS21: (1 | Gsp) + (1 | Ref_name/Ref)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS6 5 191.90 209.72 -90.949 181.90
## modSS9 5 209.77 227.59 -99.886 199.77 0.000 0 1
## modSS12 6 193.08 214.47 -90.540 181.08 18.692 1 1.536e-05 ***
## modSS15 6 206.02 227.41 -97.011 194.02 0.000 0 1
## modSS18 6 186.73 208.12 -87.364 174.73 19.294 0 < 2.2e-16 ***
## modSS21 7 188.73 213.68 -87.364 174.73 0.000 1 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: anymortSS2
## Models:
## modSS16: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 |
## modSS16: Ref_name)
## modSS17: Binary_effect_mortality ~ Sed_level_stand_mg + Sed_exposure_days +
## modSS17: (1 | Gsp) + (1 | Ref_name)
## modSS18: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days +
## modSS18: (1 | Gsp) + (1 | Ref_name)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modSS16 4 216.18 230.44 -104.089 208.18
## modSS17 5 195.39 213.21 -92.695 185.39 22.788 1 1.809e-06 ***
## modSS18 6 186.73 208.12 -87.364 174.73 10.661 1 0.001094 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Binary_effect_mortality ~ log10_conc * Sed_exposure_days + (1 |
## Gsp) + (1 | Ref_name)
## Data: anymortSS2
##
## AIC BIC logLik deviance df.resid
## 185.1 206.5 -86.6 173.1 255
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6218 -0.3392 -0.1546 -0.0635 4.9901
##
## Random effects:
## Groups Name Variance Std.Dev.
## Gsp (Intercept) 4.036 2.009
## Ref_name (Intercept) 5.147 2.269
## Number of obs: 261, groups: Gsp, 21; Ref_name, 11
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.58858 1.87960 -3.505 0.000456 ***
## log10_conc 0.74576 0.67171 1.110 0.266894
## Sed_exposure_days 0.02309 0.02624 0.880 0.378980
## log10_conc:Sed_exposure_days 0.02746 0.01556 1.765 0.077639 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Effect estimates
## [1] 0.8773946##
## p 0 1
## 0 212 25
## 1 7 17## Setting levels: control = 0, case = 1## Setting direction: controls < cases## Area under the curve: 0.9086## R2m R2c
## theoretical 0.1683615 0.7806518
## delta 0.1302107 0.6037558
## $Gsp
##
## $Ref_name
## Groups Name Std.Dev.
## Gsp (Intercept) 2.0091
## Ref_name (Intercept) 2.2687## [1] 9.666769
There is no significant relationship between exposure concentration of suspended sediment and the odds of any tissue mortality (GLMM z = 1.110, p = 0.267), nor between exposure duration and the odds of any tissue mortality (GLMM z = 0.880, p = 0.379) after accounting for variability among articles and species.
Prediction figures
#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb18_log <- predict(modSSb18_log, newdata=anymortSS2, type="response")
anymortSS3 <- cbind(anymortSS2, pred_modSSb18_log)
anymortSS_plot <- ggplot(data = anymortSS3) +
geom_point(mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_mortality,
color = Ref)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Bleaching due to sediment exposure",
color = "Study") +
scale_x_log10(limits=c(0.1,max(anymortSS3$Sed_level_stand_mg)),breaks=c(0.01,0.1,1,10,100,1000),
label=c("0.01","0.1","1","10","100","1000")) +
geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb18_log), inherit.aes=FALSE) +
theme_bw()
anymortSS_plot
That plot is a bit…confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/
By sediment exposure concentration
jvalues <- with(anymortSS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
anymortSS2$log10_conc <- j
predict(modSSb18_log, newdata = anymortSS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat4_2 <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in log10_conc values and convert to data frame
plotdat4_2 <- as.data.frame(cbind(plotdat4_2, jvalues))
plotdat4_2 <- plotdat4_2 %>% mutate(Sed_level_linear=10^jvalues)
# better names and show the first few rows
colnames(plotdat4_2) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "log10_mg_L", "mg_L")
#head(plotdat4_2)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat4_2$mg_L, plotdat4_2$PredictedProbability, xout=x))
print(round(unlist(sed2), digits=3))## x y x y x y x y
## 1.000 0.019 5.000 0.055 10.000 0.082 50.000 0.170
## x y x y x y
## 100.000 0.217 500.000 0.323 1000.000 0.369# plot average marginal predicted probabilities
ggplot(plotdat4_2, aes(x = log10_mg_L)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat4_2, aes(x = log10_mg_L)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat4_2, aes(x = mg_L, y = PredictedProbability)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1)) +
scale_x_log10()
#Overlaying average marginal predicted probabilities on figure
#by Ref
anymortSS_plot2 <- ggplot() +
geom_point(data = anymortSS2, mapping = aes(
x = Sed_level_stand_mg,
y = Binary_effect_mortality)) +
labs(x = "Sediment exposure concentration (mg/L)",
y = "Predicted probability of any tissue\nmortality to sediment exposure",
linetype = "Predicted\nProbability") +
scale_x_log10(limits=c(0.1,max(anymortSS2$Sed_level_stand_mg)),
breaks=c(0.01,0.1,1,10,100,1000,loaelSS2),
label=c("0.01","0.1","1","10","100","1000",round(loaelSS2,digits=1))) +
geom_ribbon(data = plotdat4_2, aes(x = mg_L, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat4_2, aes(x = mg_L, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat4_2, aes(x = mg_L, y = MedianProbability,
linetype = "solid")) +
geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
anymortSS_plot2
By sediment exposure duration
jvalues <- with(anymortSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))
# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
anymortSS2$Sed_exposure_days <- j
predict(modSSb18_log, newdata = anymortSS2, type = "response")
})
# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat4_2b <- t(sapply(pp, function(x) {
c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))
# add in Sed_exposure_days values and convert to data frame
plotdat4_2b <- as.data.frame(cbind(plotdat4_2b, jvalues))
# better names and show the first few rows
colnames(plotdat4_2b) <- c("PredictedProbability", "MedianProbability",
"LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat4_2b)
# plot average marginal predicted probabilities
ggplot(plotdat4_2b, aes(x = Sed_dur)) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
ggplot(plotdat4_2b, aes(x = Sed_dur)) +
geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
geom_line(aes(y = PredictedProbability), size = 2) +
geom_line(aes(y = MedianProbability), size = 0.5) +
ylim(c(0, 1))
#Overlaying average marginal predicted probabilities on figure
#by Ref
anymortSS_plot3 <- ggplot() +
geom_point(data = anymortSS2, mapping = aes(
x = Sed_exposure_days,
y = Binary_effect_mortality,
color = Ref)) +
labs(x = "Sediment exposure duration (days)",
y = "Predicted probability of any tissue\nmortality to sediment exposure",
color = "Binary Data\nby Study",
linetype = "Predicted\nProbability") +
scale_x_continuous(limits=c(0,max(anymortSS2$Sed_exposure_days))) +
geom_ribbon(data = plotdat4_2b, aes(x = Sed_dur, y = PredictedProbability,
ymin = LowerQuantile, ymax = UpperQuantile),
alpha = 0.15) +
geom_line(data = plotdat4_2b, aes(x = Sed_dur, y = PredictedProbability,
linetype = "twodash")) +
geom_line(data = plotdat4_2b, aes(x = Sed_dur, y = MedianProbability,
linetype = "solid")) +
theme_bw() +
scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))
anymortSS_plot3